Image Deconvolution using Sparse Regularization

This numerical tour explores the use of sparsity regularization to perform image deconvolution.

Contents

Installing toolboxes and setting up the path.

You need to download the following files: signal toolbox and general toolbox.

You need to unzip these toolboxes in your working directory, so that you have toolbox_signal and toolbox_general in your directory.

For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.

Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab command you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.

Execute this line only if you are using Matlab.

getd = @(p)path(p,path); % scilab users must *not* execute this

Then you can add the toolboxes to the path.

getd('toolbox_signal/');
getd('toolbox_general/');

Sparse Regularization

This tour consider measurements y=Φf0+w$y=\Phi f_0 + w$ where Φ$\Phi$ is a convolution Φf=h⋆h$\Phi f = h \star h$ and w$w$ is an additive noise.

This tour is focussed on using sparsity to recover an image from the measurements y$y$. It consider a synthesis-based regularization, that compute a sparse set of coefficients (a⋆m)m$(a_m^{\star})_m$ in a frame Ψ=(ψm)m$\Psi = (\psi_m)_m$ that solves

$$a^{\star} \in \text{argmin}_a \: \frac{1}{2}\|y-\Phi \Psi a\|^2 + \lambda J(a)$$

where λ$\lambda$ should be adapted to the noise level ∥w∥$\|w\|$

Here we used the notation

$$\Psi a = \sum_m a_m \psi_m$$ to indicate the reconstruction operator, and J(a)$J(a)$ is the ℓ1$\ell^1$ sparsity prior $$J(a)=\sum_m \|a_m\|.$$

Image Blurring

Deconvolution corresponds to removing a blur from an image. We use here a Gaussian blur.

Parameters for the tour: width of the kernel (in pixel) and noise level.

setting = 1;
switch setting
    case 1
        % difficult
        s = 3;
        sigma = .02;
    case 2
        % easy
        s = 1.2;
        sigma = .02;
end

First we load the image to be processed.

n = 128*2;
name = 'lena';
name = 'boat';
name = 'mri';
f0 = load_image(name);
f0 = rescale(crop(f0,n));

Display it.

clf;
imageplot(f0);

We build a convolution kernel. Since we are going to use Fourier to compute the convolution, we set the center of the kernel in the (1,1) pixel location.

Kernel.

x = [0:n/2-1, -n/2:-1];
[Y,X] = meshgrid(x,x);
h = exp( (-X.^2-Y.^2)/(2*s^2) );
h = h/sum(h(:));

Useful for later : the Fourier transform (should be real because of symmetry).

hF = real(fft2(h));

Display the kernel h$h$ and its transform hˆ$\hat h$. We use fftshift to center the filter for display.

clf;
imageplot(fftshift(h), 'Filter', 1,2,1);
imageplot(fftshift(hF), 'Fourier transform', 1,2,2);

We use this short hand for the filtering. Scilab user should define a function in a separate file to perform this. Note that this is a symmetric operator.

Phi = @(x)real(ifft2(fft2(x).*hF));

Apply the filter.

y0 = Phi(f0);

Display the filtered observation.

clf;
imageplot(f0, 'Image f0', 1,2,1);
imageplot(y0, 'Observation without noise', 1,2,2);

Add some noise of variance σ2$\sigma^2$, to obtain y=Φf0+w=f0⋆h+w$y=\Phi f_0 + w = f_0 \star h + w$.

y = y0 + randn(n,n)*sigma;

Display.

clf;
imageplot(y0, 'Observation without noise', 1,2,1);
imageplot(clamp(y), 'Observation with noise', 1,2,2);

Soft Thresholding in a Basis

The soft thresholding operator is at the heart of ℓ1$\ell^1$ minimization schemes. It can be applied to coefficients a$a$, or to an image f$f$ in an ortho-basis.

The soft thresholding is a 1-D functional that shrinks the value of coefficients.

$$s_T(u)=\max(0,1-T/|u|)u$$

Define a shortcut for this soft thresholding 1-D functional.

SoftThresh = @(x,T)x.*max( 0, 1-T./max(abs(x),1e-10) );

Display a curve of the 1D soft thresholding.

clf;
T = linspace(-1,1,1000);
plot( T, SoftThresh(T,.5) );
axis('equal');

Note that the function SoftThresh can also be applied to vector (because of Matlab/Scilab vectorialized computation), which defines an operator on coefficients:

$$S_T(a) = ( s_T(a_m) )_m.$$

In the next section, we use an orthogonal wavelet basis Ψ$\Psi$.

We set the parameters of the wavelet transform.

Jmax = log2(n)-1;
Jmin = Jmax-3;

Shortcut for Ψ$\Psi$ and Ψ∗$\Psi^*$ in the orthogonal case.

options.ti = 0; % use orthogonality.
Psi = @(a)perform_wavelet_transf(a, Jmin, -1,options);
PsiS = @(f)perform_wavelet_transf(f, Jmin, +1,options);

The soft thresholding opterator in the basis Ψ$\Psi$ is defined as

$$S_T^\Psi(f) = \sum_m s_T( \langle f,\psi_m \rangle ) \psi_m$$

It thus corresponds to applying the transform Ψ∗$\Psi^*$, thresholding the coefficients using ST$S_T$ and then undoing the transform using Ψ$\Psi$.

$$S_T^\Psi(f) = \Psi \circ S_T \circ \Psi^*$$

SoftThreshPsi = @(f,T)Psi(SoftThresh(PsiS(f),T));

This soft thresholding corresponds to a denoising operator.

clf;
imageplot( clamp(SoftThreshPsi(f0,.1)) );

Deconvolution using Orthogonal Wavelet Sparsity

If Ψ$\Psi$ is an orthogonal basis, a change of variable shows that the synthesis prior is also an analysis prior, that reads

$$f^{\star} \in \text{argmin}_f \: E(f) = \frac{1}{2}\|y-\Phi f\|^2 + \lambda \sum_m \|\langle f,\psi_m \rangle\|.$$

To solve this non-smooth optimization problem, one can use forward-backward splitting, also known as iterative soft thresholding.

It computes a series of images f(ℓ)$f^{(\ell)}$ defined as

$$f^{(\ell+1)} = S_{\tau\lambda}^{\Psi}( f^{(\ell)} - \tau \Phi^{*} (\Phi f^{(\ell)} - y) )$$

Set up the value of the threshold.

lambda = .02;

In our setting, since h$h$ is symmetric, one has Φ∗f=Φf=f⋆h$\Phi^* f = \Phi f = f \star h$.

For f(ℓ)$f^{(\ell)}$ to converge to a solution of the problem, the gradient step size should be chosen as

$$\tau < \frac{2}{\|\Phi^* \Phi\|}$$

Since the filtering is an operator of norm 1, this must be smaller than 2.

tau = 1.5;

Number of iterations.

niter = 100;

Initialize the solution.

fSpars = y;

First step: perform one step of gradient descent of the energy ∥y−f⋆h∥2$\|y-f\star h\|^2$.

fSpars = fSpars + tau * Phi( y-Phi(fSpars) );

Second step: denoise the solution by thresholding.

fSpars = SoftThreshPsi( fSpars, lambda*tau );

Exercice 1: (the solution is exo1.m) Perform the iterative soft thresholding. Monitor the decay of the energy E$E$ you are minimizing.

exo1;

Display the result.

clf;
imageplot(clamp(fSpars), ['Sparsity deconvolution, SNR=' num2str(snr(f0,fSpars),3) 'dB']);

Exercice 2: (the solution is exo2.m) Try to find the best threshold λ$\lambda$. To this end, perform a lot of iterations, and progressively decay the threshold λ$\lambda$ during the iterations. Record the best result in fBestOrtho.

exo2;

Display the result.

clf;
imageplot(clamp(fBestOrtho), ['Sparsity deconvolution, SNR=' num2str(snr(f0,fBestOrtho),3) 'dB']);

Deconvolution using Translation Invariant Wavelet Sparsity

Orthogonal sparsity performs a poor regularization because of the lack of translation invariance. This regularization is enhanced by considering Ψ$\Psi$ as a redundant tight frame of translation invariant wavelets.

One thus looks for optimal coefficients a⋆$a^\star$ that solves

$$a^{\star} \in \text{argmin}_a \: \frac{1}{2}\|y-\Phi \Psi a\|^2 + \lambda J(a)$$

One should be careful that because of the redundancy of the wavelet tight frame, one should use a weighted ℓ1$\ell^1$ norm, where each coefficient is divided by the number of redundancy at each scale.

$$J(a) = \sum_{j,k} u_j \|a_{j,k}\|$$

where the wavelet coefficients am=aj,k$a_{m} = a_{j,k}$ are indexed by the scale (and orientation) j$j$ and the location k$k$.

Compute the scaling factor (inverse of the redundancy).

J = Jmax-Jmin+1;
u = [4^(-J) 4.^(-floor(J+2/3:-1/3:1)) ];

Value of the regularization parameter.

lambda = .01;

Shortcut for the wavelet transform. Important: Note that PsiS is the shortcut for Ψ∗$\Psi^*$, but Psi is the shortcut for Ξ=(Ψ∗)+$\Xi = (\Psi^*)^+$ that satisfy ΞΦ∗f=f$\Xi \Phi^* f = f$.

options.ti = 1; % use translation invariance
Psi = @(a)perform_wavelet_transf(a, Jmin, -1,options);
PsiS = @(f)perform_wavelet_transf(f, Jmin, +1,options);

The forward-backward algorithm now compute a series of wavelet coefficients a(ℓ)$a^{(\ell)}$ computed as

$$a^{(\ell+1)} = S_{\tau\lambda}( a^{(\ell)} + \Psi^*\Phi( y - \Phi\Xi a^{(\ell)} ) ).$$

The soft thresholding is defined as:

$$\forall m, \quad S_T(a)_m = \max(0, 1-T/\|a_m\|)a_m.$$

The step size should satisfy:

$$\tau < \frac{2}{\|\Phi\Psi\|^2}.$$

tau = 1.5;

Initialize the wavelet coefficients with those of the observations.

a = PsiS(y);

Gradient descent.

a = a + tau * PsiS( Phi( y-Phi(Psi(a)) ) );

Soft threshold.

a = SoftThresh( a, lambda*tau );

Important: keep in mind that the prior J(a)$J(a)$ is a weighted ℓ1$\ell^1$ norm, it should thus be computed this way:

U = repmat( reshape(u,[1 1 length(u)]), [n n 1] );
Ja = sum(sum(sum( abs(a.*U) )));

Exercice 3: (the solution is exo3.m) Perform the iterative soft thresholding. Monitor the decay of the energy.

exo3;

Perform the reconstruction.

fTI = Psi(a);

Display the result.

clf;
imageplot(fTI);

Exercice 4: (the solution is exo4.m) Compute the optimal value of λ$\lambda$, and record the optimal reconstruction fBestTI.

exo4;

Display the result.

clf;
imageplot(clamp(fBestTI), ['Sparsity deconvolution TI, SNR=' num2str(snr(f0,fBestTI),3) 'dB']);

Exercice 5: (the solution is exo5.m) Compare with the result of TV regularization, record the optimal TV result in fBestTV.

exo5;

Display the result.

clf;
imageplot(clamp(fBestTV), ['TV deconvolution, SNR=' num2str(snr(f0,fBestTV),3) 'dB']);

Copyright (c) 2010 Gabriel Peyre